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Abstract 

As sensors and flow control actuators become smaller, cheaper, and more pervasive, the 
use of feedback control to manipulate the details of fluid flows becomes increasingly attractive. 
One of the challenges is to develop mathematical models that describe the fluid physics relevant 
to the task at hand, while neglecting irrelevant details of the flow in order to remain computa- 
tionally tractable. A number of techniques are presently used to develop such reduced-order 
models, such as proper orthogonal decomposition (POD), and approximate snapshot-based bal- 
anced truncation, also known as balanced POD. Each method has its strengths and weaknesses: 
for instance, POD models can behave unpredictably and perform poorly, but they can be com- 
puted directly from experimental data; approximate balanced truncation often produces vastly 
superior models to POD, but requires data from adjoint simulations, and thus cannot be applied 
to experimental data. 

In this paper, we show that using the Eigensystem Realization Algorithm (ERA) lfl5ll . one 
can theoretically obtain exactly the same reduced order models as by balanced POD. Moreover, 
the models can be obtained directly from experimental data, without the use of adjoint infor- 
mation. The algorithm can also substantially improve computational efficiency when forming 
reduced-order models from simulation data. If adjoint information is available, then balanced 
POD has some advantages over ERA: for instance, it produces modes that are useful for mul- 
tiple purposes, and the method has been generalized to unstable systems. We also present a 
modified ERA procedure that produces modes without adjoint information, but for this pro- 
cedure, the resulting models are not balanced, and do not perform as well in examples. We 
present a detailed comparison of the methods, and illustrate them on an example of the flow 
past an inclined flat plate at a low Reynolds number. 
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1 Introduction 



In the last decade, substantial developments have been made in the area of model-based feedback 
flow control of fluids: for instance, see the recent reviews by ifTTl 171181. In many applications, the 
focus is on how to apply actuation to maintain the flow around an equilibrium state of interest, for 
instance to delay transition to turbulence, or control separation on a bluff body. Linear control theory 
provides efficient tools for the analysis and design of feedback controllers. However, a significant 
challenge is that models for flow control problems are often very high dimensional, e.g., on the order 
of O(10 5 ~ 9 ), so large that it becomes computationally infeasible to apply linear control techniques. 
To address this issue, model reduction, by which a low-order approximate model is obtained, is 
therefore widely employed. 

Several techniques are available for model reduction, many of which involve projection onto a set 
of modes. These may be global eigenmodes of a linearized operator [1], modes determined by 
proper orthogonal decomposition (POD) of a set of data [,' 131 , and various variants of POD, such 
as including shift modes ETTl . An approach that is used widely for model reduction of linear sys- 
tems is balanced truncation [20], and while this method is computationally intractable for systems 
with very large state spaces (dimension > 10 5 ), recently an algorithm for computing approximate 
balanced truncation from snapshots of linearized and adjoint simulations has been developed l23l 
and successfully applied to a variety of high-dimensional flow control problems Ifl4l l2ll4ll. In this 
method, sometimes called balanced POD, one obtains two sets of modes (primal and adjoint) that 
are bi-orthogonal, and uses those for projection of the governing equations, just as in standard POD. 
Compared to most other methods, including POD, balanced truncation has key advantages, such as 
a priori error bounds, and guaranteed stability of the reduced-order model (if the original high-order 
system is stable). Balanced POD is an approximation of exact balanced truncation that is compu- 
tationally tractable when the number of states is very large (for instance, up to 10 7 ), and typically 
produces models that are far more accurate than standard POD models. For instance, for a linearized 
channel flow investigated in fPfll . even though the first 5 POD modes capture over 99.7% of the en- 
ergy in a dataset exhibiting large transient growth, a low-dimensional model obtained by projection 
onto these modes completely misses the transient growth. By contrast, a 3-mode balanced POD 
model captures the transient growth nearly perfectly; to do as well with a standard POD model, 17 
modes were required. 

The main steps of balanced POD include (a) taking snapshots from impulse responses of the lin- 
earized and adjoint systems, (b) computing a singular value decomposition (SVD) of a matrix 
formed from inner products of these snapshots, (c) constructing primal modes and adjoint modes 
from the resulting singular vectors, and (d) projecting the high-dimensional dynamics onto these 
modes. 

While effective in many examples, balanced POD also faces challenges, especially for use with ex- 
perimental data. The main restriction is that balanced POD requires snapshots of impulse-response 
data from an adjoint system, and adjoint information is not available for experiments. 
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To address this issue, here we describe an algorithm widely used for system identification and 
model reduction, the eigensystem realization algorithm (ERA) IPT51 . This algorithm has been used 
for problems in fluid mechanics, primarily as a system-identification technique for flow control flU 
13, but also for model reduction iTTTl l24ll . Our main result, presented in Section [2| is that, for 
linear systems, ERA theoretically produces exactly the same reduced-order models as balanced 
POD, with no need of an adjoint system, and at an order of magnitude lower computational cost. 
This result implies that one can realize approximate balanced truncation even in experiments, and 
can also improve computational efficiency in simulations. We note that ERA and snapshot-based 
approximate balanced truncation have been applied together in a model reduction procedure in ifTOl . 
However, the theoretical equivalence between these two algorithms was not explored in that work. 

We present a comparison between balanced POD and ERA, and show that if adjoint information 
is available, balanced POD also has its own advantages. In particular, balanced POD provides 
sets of bi-orthogonal primal/adjoint modes for the linear system, and can be directly generalized 
to unstable systems. In Section [3] we discuss a modified ERA algorithm that, in the absence of 
adjoint simulations, uses "pseudo- adjoint modes" to compute reduced order models; however, this 
method does not produce balanced models, and performs worse than balanced POD in examples. In 
Section |4| we illustrate these methods using a numerical example of the two-dimensional flow past 
an inclined plate, at a low Reynolds number. 

2 The eigensystem realization algorithm as snapshot-based approxi- 
mate balanced truncation 

In this section, we summarize the steps involved in approximate balanced truncation (balanced 
POD), and the Eigensystem Realization Algorithm, and show that they are equivalent. 

Balanced truncation involves first constructing a a coordinate transformation that "balances" a lin- 
ear input-output system, in the sense that certain measures of controllability and observability (the 
Gramian matrices) become diagonal and identical ll20ll . A reduced-order model is then obtained 
by truncating the least controllable and observable states, which correspond to the smallest diag- 
onal entries in the transformed system. Unfortunately, the exact balanced truncation algorithm is 
not tractable for the large state dimensions encountered in fluid mechanics. However, an approx- 
imate, snapshot-based balanced truncation algorithm, referred to as Balanced Proper Orthogonal 
Decomposition (balanced POD) was proposed in ||23l , and has been used successfully in several 
examples lPT4ll2ll4l. 

The second technique, the eigensystem realization algorithm (ERA), has been used both for sys- 
tem identification and for model reduction, and it is well known that the models produced by ERA 
are approximately balanced lfl2l [Toll . Here we show further that, theoretically, ERA produces ex- 
actly the same reduced order models as balanced POD. This equivalence indicates that ERA can 
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be regarded as an approximate balanced truncation method, in the sense that, before truncation, it 
implicitly realizes a coordinate transformation under which a pair of approximate controllability 
and observability Gramians are exactly balanced. This feature distinguishes ERA from other model 
reduction methods that first realize truncations and then balance the reduced order models. Note 
that in ERA the Gramians, and the balancing transformation itself, are never explicitly calculated, 
as we will also show in the following discussions. 

For both techniques, we will consider a high-dimensional, stable, discrete-time linear system, de- 
scribed by 

x(k + 1) = Ax(k) + Bu(k) 

y(k) = Cx(k), (1) 

where k G Z is the time step index, u(k) G W denotes a vector of inputs (for instance, actuators 
or disturbances), y(k) G M. q a vector of outputs (for instance, sensor measurements, or simply 
quantities that one wishes to model), and x(k) G M" denotes the state variable (for instance, flow 
variables at all gridpoints of a simulation). These equations may arise, for instance, by discretizing 
the Navier-Stokes equations in time and space, and linearizing about a steady solution, as will be 
demonstrated in the example in Section|4] The goal is to obtain an approximate model that captures 
the same relationship between inputs u and outputs y, but with a much smaller state dimension: 

x r {k + 1) = A r x r {k) + B r u{k) 
y{k) = C r x r {k) 

where the reduced state variable x r (k) G W, r -C n. We consider the discrete-time setting, because 
we are primarily interested in discrete-time data from simulations or experiments. 



2.1 Snapshot- based approximate balanced truncation (balanced POD) 

Here, we give only a brief overview of the balanced POD algorithm, and for details of the method, 
we refer the reader to (23). The algorithm involves three main steps: 



Step 1: Collect snapshots. Run impulse-response simulations of the primal system ([T]) and 
collect m c + 1 snapshots of states x(k) in m c P + 1 steps: 

X=[B A P B A 2P B ■ ■ ■ A mcP B] , (3) 

where P is the sampling period. In addition, run impulse-response simulations for the adjoint 
system 

z(k + l) = A*z(k) + C*v(k) (4) 

where the asterisk * stands for adjoint of a matrix, and collect m a + 1 snapshots of states z(k) 
in m a P + 1 steps: 



Y 



C* (A*) P C* {A*) 2P C*--- (A*) m ° p C* 



(5) 
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Calculate the generalized Hankel matrix, 

H = Y*X. 



(6) 



Step 2: Compute modes. Compute the singular value decomposition of H: 



H = UY,V* 



[Ui U 2 



"Si 0" 




~v*~ 







V* 



(7) 



where the diagonal matrix Si € I" 1 xni is invertible and includes all non-zero singular values 
of H, m = rahk(H), and U*U\ = V*V\ = I m xm- Choose r < n\. Let U r and V r denote 
the sub-matrices of U\ and V\ that include their first r columns, and S r the first r x r diagonal 
block of Si. Calculate 



XV r X r 5 ; 



(8) 



where the columns of <£ r and ^ r are respectively the first r primal and adjoint modes of 
system ([T]). The two sets of modes are bi-orthogonal: ^*$ r = / rxr . 

Step 3: Project dynamics. The system matrices in the reduced order model ([2]) are 



A r = V*A<l> r ; B r = V r B; C r = C$ r . 



(9) 



Note that the n x n controllability/observability Gramians are approximated by the matrices XX* 
and YY*. The reduced-order model ([2]) is obtained by considering a subspace x = $ r x r , and 
projecting the dynamics ([T]) onto this subspace using the adjoint modes given by *I' r . It was shown 
in If23! that <£ r and ^ r respectively form the first r columns of the balancing transformation/inverse 
transformation that exactly balance the approximate controllability/observability Gramians XX* 
and YY*; see more discussion in Section[3] 



2.2 The eigensystem realization algorithm 



The eigensystem realization algorithm (ERA) was proposed in lfl3Tl as a system identification and 
model reduction technique for linear systems. The algorithm follows three main steps |[T5l[T6ll : 



• Step 1: Run impulse -response simulations/experiments of the system ([!]) for {m c + m )P + 2 
steps, where m c and m respectively reflect how much effect is taken for considering con- 
trollability and observability, and P again is the sampling period. Collect the snapshots of the 
outputs y in the following pattern: 

(CB, CAB, CA P B, CA P+1 B, ... 

(JA m c p B CA mcP+l B (jj^(m c +m )P ^(m t+mo )P+l^j 
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The terms CA k B are commonly called Markov parameters. Construct a generalized Hankel 
matrix H £ R9K,+l)xp(m c +l) 



H 



CB 
CA P B 



CA P B 
CA 2P B 



CA m ° p B CA( m °+V p B 



CA m - p B 

CA (m c +l)P_g 
Q jS L (m c +m )P £> 



(ID 



Step 2: Compute SVD of H, exactly as in Q, to obtain U l3 V\, Ej. Let r < rank(if). Let 
C/ r and y r denote the sub-matrices of U\ and Vi that include their first r columns, and E r the 
first r x r diagonal block of Ei. 

Step 3: The reduced A r , B T and C r in ([2]) are then defined as 



Arf — Ey 2 t/^. -/^ V^-E^ 2 \ 

i 

Z? r = the first p columns of E 2 V* 

i 

C r = the first q rows of L^E 2 



where 



H' 



CAB 



CA p+l B 



CA mcP+1 B 



(12) 



(13) 



which can again be constructed directly from the collected snapshots ( 10 1 



2.3 Theoretical equivalence between ERA and balanced POD 



The first observation is that, with X and Y given by Q and ([5]), the generalized Hankel matrices 
obtained in balanced POD and ERA, respectively by Q and ( [TTj ), are theoretically identical. The 
theoretical equivalence between the two algorithms then follows immediately: First, H' given in 
(13 1 satisfies H' = Y*AX, which implies the matrices A r obtained in the two algorithms are 

identical. To show the equivalence of B r , first note that the SVD hi leads to T,^ 2 U*H = E^Vf, 

_ i i 

which, by definition of U r , V r , E r , implies E r 2 U*H = HfiV*. (Note that it does not imply 

H = U r Z r V r *, since U r U* is not the identity.) Thus, in balanced POD, B r = %B = E^U*Y*B, 

_i i 

which equals the first p columns of E r 2 U*H = Er V* , which is the B r given by ERA. Similarly, 

the SVD 111 leads to HViE^ 5 = ETjE? and then HV r Y>r* = U r 2$. Thus, in balanced POD, 

_I _i I 

C r = C&r = CXV r T, r 2 , which equals the first q rows of HV r Yl r 2 = U r T, 2 , the C r given by 

ERA. In summary, we have: 
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Main result. The reduced system matrices A r , B r and C r generated in balanced POD and ERA, 
respectively by ([9]) and (12 \, are theoretically identical. 



In practice, these two algorithms may generate slightly different reduced order models, because the 
Hankel matrices calculated in the two algorithms are usually not exactly the same, due to small nu- 
merical inaccuracies in adjoint simulations, and/or in matrix multiplications needed to compute the 
sub-blocks in the Hankel matrices. In the following discussions, we compare these two algorithms 
in more detail. 



2.4 Comparison between ERA and balanced POD 



While ERA and balanced POD produce theoretically identical reduced-order models, the techniques 
differ in several important ways, both conceptually and computationally. Neither ERA nor balanced 
POD calculate Gramians explicitly, but balanced POD does construct approximate controllability 
and observability matrices X and Y* , from which one calculates the generalized Hankel matrix H 
and balancing transformation. Balanced POD thus incurs additional computational cost, because 
one needs to construct the adjoint system ([4]), run adjoint simulations for Y, and then calculate each 
block of H by matrix multiplication. Thus we see that the advantages of ERA include: 



1. Adjoint-free: ERA is a feasible balanced truncation method for experiments, since it needs 
only the output measurements from the response to an impulsive input. Note that ERA has 
been successfully applied in several flow control experiments [ 6 , 5 ] , as a system- identification 
technique rather than a balanced-truncation method. In practice, input-output sensor re- 
sponses are often collected by applying a broadband signal to the inputs, and the ARMARKOV 
method |3][TH can then be used to identify the Markov parameters, or even directly the gen- 
eralized Hankel matrix, from the input-output data history. 

2. Computational efficiency: For large problems, typically the most computationally expen- 
sive component of computing balanced POD is constructing the generalized Hankel matrix H 
in ([6]), as this involves computing inner products of all of the (large) primal and adjoint snap- 



shots with each other. ERA is significantly more efficient at constructing the matrix H in ( 1 1 
since only the first row and last column of block matrices, i.e., the (m c + m + 1) Markov 
parameters, need be obtained by matrix multiplication. All the other m c x m block matrices 
in H are copies of other blocks, and need not be recomputed. For balanced POD, the matrix H 
is obtained by computing all the (m c + 1) x (m + 1) matrix multiplications (inner products) 
between corresponding blocks in Y* and X in ([6]>. Thus, for example, if m c = tuq = 200, the 
computing time needed for constructing H in ERA will be about only 1% of that in balanced 
POD. See Table [T] for a detailed comparison on computational efficiency between balanced 
POD and ERA in the example of the flow past an inclined flat plate. 



At the same time, balanced POD also provides its own advantages: 
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1. Sets of bi-orthogonal primal/adjoint modes: Balanced POD provides sets of bi-orthogonal 
primal/adjoint modes, the columns of $ r and In comparison, without the adjoint sys- 
tem, ERA cannot provide the primal and adjoint modes. At best, the primal modes may be 
computed, using the first equation in ([8]), if the matrix X Q is stored (in addition to the 
Markov parameters). But the adjoint modes cannot be computed without solutions of the 
adjoint system. In this sense, balanced POD incorporates more of the physics of the system 
(the two sets of bi-orthogonal modes), while ERA is purely based on input-output data of 
the system. The primal/adjoint modes together can be useful for system analysis and con- 
troller/observer design purposes in several ways: for instance, in flow control applications, a 
large-amplitude region from the most observable mode (the leading adjoint mode) can be a 
good location for actuator placement. Also, although balanced POD is a linear method, a non- 
linear system can be projected onto these sets of modes to obtain a nonlinear low-dimensional 
model. For instance, the transformation x = $ r x r , x r = ty*x can be employed to reduce a 
full-dimensional nonlinear model x = f(x) to a low-dimensional system x r = **/($ r » r ). 
Finally, if parameters (such as Reynolds number or Mach number) are present in the original 
equations, balanced POD can retain these parameters in the reduced-order models. When the 
values of parameters change, the reduced order model by balanced POD may still be valid 
and perform well; see lPl4l for an application to linearized channel flow. 

2. Unstable systems: Balanced POD has been extended to neutrally stable [ 19] and unstable 
systems 0. In those cases, one first calculates the right/left eigenvectors corresponding to 
the neutral/unstable eigenvalues of the state-transition matrix A, using direct/adjoint simu- 
lations. Using these eigenvectors, the system is projected onto a stable subspace and then 
balanced truncation is realized for the stable subsystem. ERA for unstable systems is still an 
open problem, if adjoint operators are not available. However, we note that, once the stable 
subsystem is obtained, ERA can still be applied to it and efficiently realize its approximate 
balanced truncation. 



ERA for systems with high-dimensional outputs. The method of output projection proposed in 
l23l makes it computationally feasible to realize approximate balanced truncation for systems with 
high-dimensional outputs — for instance, if one wishes to model the entire state x, say the flow field 
in the entire computational or experimental domain. This method involves projecting the outputs 
onto a small number of POD modes, determined from snapshots of y from the impulse-response 
dataset. This method can be directly incorporated into ERA as follows: First, run impulse response 
simulations of the original system and collect Markov parameters as usual. Then, compute the 
leading POD modes of the dataset of Markov parameters and stack them as columns of a matrix 0. 
Left multiply those Markov parameters by 0* to project the outputs onto these POD modes. A 
generalized Hankel matrix is then constructed using these modified Markov parameters, and the 
usual steps of ERA follow. 
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3 A modified ERA method using pseudo-adjoint modes 



We have seen that one of the drawbacks of ERA is that it does not provide modes that could be 
used, for instance, for projection of nonlinear dynamics, or to retain parameters in the models. 

— 1/2 

More precisely, using ERA, one may still obtain primal modes $1 = XV\T, 1 as in balanced 
POD (see (|7]-[8]>), as long as the state snapshots are collected and stored in X. But it is not possible 
to obtain the corresponding adjoint modes necessary for projection, without performing adjoint 
simulations to gather snapshots for the matrix Y. This is a severe drawback, as adjoint solutions can 
be expensive to perform, and are not available for experimental data. One idea, proposed in 11221 . is 
to define a set of approximate adjoint modes using the Moore-Penrose pseudo-inverse of <&i: 

*i = SiC*?*!) -1 . (14) 

We will call the adjoint modes as defined above the pseudo-adjoint modes corresponding to the 
modes $x- The system matrices of a r-dimensional reduced-order model (r < rank(i^)) generated 
by this approach then read 

A r = %A$ r - B r = y* r B; C r = C<Z> r , (15) 
where <£ r and are respectively the first n x r sub-blocks of $1,^1. 

While this idea does produce a set of modes that can be used for projection, we show in this sec- 
tion that, unfortunately, the resulting transformation is not a balancing transformation, and does not 
produce models that are an approximation to balanced truncation. In fact, the resulting models are 
closer to those produced by the the standard POD/Galerkin method: as with standard POD/Galerkin, 
the method performs well as long as the most controllable and most observable directions coincide. 
However, when these directions differ (as is the case for many problems of interest, including the 
example in Section the method performs poorly. These systems in which controllable and ob- 
servable directions do not coincide are precisely the systems for which balanced POD and ERA give 
the most improvement over the more traditional POD/Galerkin approach. 



3.1 Transformed approximate Gramians 

First, let us recall in what sense the model -reduction procedures described in Section [2] are approx- 
imations to balanced truncation. Suppose that we have an approximation of the controllability and 
observability Gramians, factored as 

W C = XX*, W Q = YY\ (16) 

where X and Y are the matrices of snapshots from ([3]) and ([5]). In balanced POD, we define the 
primal modes as columns of $1 = XVi'E 1 ~ 2 , and the adjoint modes as columns of ^\ = YU\E 1 5 , 
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where U\, V\, and Si are defined in ([7]). We will assume in this section that the number of columns 
of X and Y (the number of snapshots, m c and m , respectively) is smaller than the number of rows 
(the state dimension, n), which is always true for the large fluids systems of interest here. 

Then balanced POD is an approximation to balanced truncation in the following sense: as shown 
in the appendix of |f23] (the proof of Proposition 2), we may construct a full (invertible, n x n) 
transformation 

r=[*i $ 2 ] (17) 

by choosing $2 such that ^*$2 = 0. That is, columns of $2 are orthogonal to the adjoint modes, 
which are columns of ^1. The inverse transformation then has the form 



(18) 



where is the matrix of adjoint modes, and ^2 is defined by (fl~8]). Then, Proposition 2 of 



states that the transformed approximate Gramians ( |T6| ) have the form 

T^WciT- 1 )* 



Si 
Mi 



T*W Q T 



Si 
M 2 



and furthermore the product of the approximate Gramians, in the transformed coordinates, is 



T^WrWnT 











(19) 



(20) 



In this sense, the transformation T balances the approximate Gramians as closely as possible: the 
Gramians are block diagonal, and the upper-left blocks are equal and diagonal. Furthermore, all 
of the states in the lower-right block (i.e., involving M\ and M2 above) are either unobservable or 
uncontrollable, as they do not appear in the product of the Gramians. 

However, if the pseudo-adjoint modes ^1 are used in place of the true adjoint modes then this 
result does not hold, as we now show. Note that, in order for the first block of rows of T- 1 to equal 
we must now define 

f = [$! $ 2 ] (21) 

where = 0. Since the range of ^1 equals the range of $1, this is then equivalent to choosing 

<l>2 such that its columns are orthogonal to the columns of $1 (the primal modes), while when the 
"true" adjoint modes are used, columns of $2 are chosen to be orthogonal to the adjoint modes 



Defining ^2 by 



*2. 



(22) 
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(b) Transformed W c , W Q and W C W Q by projection using pseudo-adjoint modes 



Figure 1: Transformed Gramian matrices: (a) using true adjoint modes (eq. ( 19-20)) and (b) using 
pseudo-adjoint modes (eq. (p3]>). Here, X and Y in ( [T6| ) are random matrices with n = 200 states 
50 snapshots. 



and m c = m, 



one can then show that, as long as rank(X) < rank(y)Q the transformed Gramians have the form 



with 



Si 
Mi 



-1/2 



T*W„T 



Si M 3 

M 3 * M 2 



T-'WcWoT 



M 3 = Si**$ 2 , 



SiM 3 
MiM| 
(23) 



(24) 



where ^1 = YU\Ti 1 are the true adjoint modes. Note that, when the true adjoint modes are used 



to define the inverse ( 18 1, then M3 = 0, since ^ / * < I ) 2 = 0. However, when pseudo-adjoint modes 
are used, then M3 is no longer zero, and in fact, can be quite large. 

An example is shown in Figure [T] which shows the magnitude of the elements of the transformed 
Gramians, where X and Y in ( [T6] > are chosen at random. Note that when true adjoint modes are 
used, the transformed Gramians are equal and diagonal, while when the pseudo- adjoint modes are 
used, the off-diagonal blocks of the transformed observability Gramian, and the product of the 
Gramians, have significant magnitude. 

'if rank(X) > rank(F), then the situation is worse, and the transformed controllability Gramian is not block 
diagonal, nor does its upper-left block equal £1. 
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Thus, when pseudo-adjoint modes are used, the resulting transformation is not, in general, a balanc- 
ing transformation: even though the upper-left blocks of the transformed Gramians are still equal 
and diagonal, the transformed observability Gramian is not block diagonal, and so its eigenvalues 
and eigenvectors do not correspond to those of the transformed controllability Gramian. Note that 
this is the whole point of balanced truncation: to transform to coordinates in which the most control- 
lable directions (dominant eigenvectors of W c ) correspond to the most observable directions (dom- 
inant eigenvectors of W ). Therefore, while the approximate balanced truncation procedure de- 
scribed in Section |2~T] exactly balances the approximate Gramians, transforming by pseudo-adjoint 
modes does not represent balancing in any meaningful sense. 

Note that the matrix M3 describes the degree to which projection using pseudo-adjoint modes fails 
to balance the approximate Gramians. This matrix equals zero if the adjoint modes (columns of 
are spanned by the primal modes (columns of $1). However, M3 is the largest when the dominant 
adjoint modes (columns of ^1) are nearly orthogonal to the dominant primal modes (columns of 
<&i). Unfortunately, this is the case in many problems of interest, in particular those involving non- 
normality: the directions spanned by the primal modes often do not coincide with the directions 
spanned by the adjoint modes. 

In the next section, we apply this approach to the flow past a flat plate, and compare it to the methods 
described in Section [2] 



4 Example: flow past an inclined flat plate 

In this section, we illustrate the application of ERA as an approximate balanced truncation method 
using a numerical example, by obtaining reduced-order models of a large-dimensional fluid system. 
We compare the resulting models with those obtained using the balanced POD method of [23 ], ERA 
with pseudo-adjoint modes as described in Section[3j and the standard POD/Galerkin method [13]. 

4.1 Model problem and parameters 

The model problem that we consider is a two-dimensional uniform flow over a flat plate inclined at 
an angle a = 25°, at a low Reynolds number Re = 100. At these conditions, the flow asymptoti- 
cally reaches a stable steady state, the streamlines of which are plotted in Figure [2] The numerical 
method used for all computations is a fast formulation of the immersed boundary method developed 
by 13, and solves for the vorticity field at each time step. We treat farfield boundary conditions us- 
ing the multiple-grid scheme described in [9] (Section 4) with five nested grids, each with 250 x 250 
points. The finest grid covers the region [—2, 3] x [—2.5, 2.5], and the largest grid covers the re- 
gion [—32,48] x [—40,40], where lengths are non-dimensionalized by the chord of the flat plate, 
whose center is located at the origin. The time step used for all simulations is 0.01 (nondimension- 
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Figure 2: Streamlines of the stable steady state past a flat plate at a = 25° (left), and the contour- 
lines of the vorticity field obtained from an impulsive input to the actuator (right). 

alized by chord and freestream velocity). The numerical model is the same as that considered in 
where balanced POD is applied for feedback controller design to stabilize an unstable steady state 
corresponding to a high angle of attack. However, here we consider the case of a stable steady state 
(with an angle of attack at 25°), for comparison of reduced order models. 



4.2 Input and output 

The governing equations are first linearized about the stable steady state, resulting in a high-dimensional 
model in the form of equation ([TJ, where the state x consists of the discrete vorticity field at the grid 
points. See (2j for the details of the linearized (and adjoint) equations and their numerical formula- 
tions. The system input u is a disturbance (or actuator) shown in Figure [2| modeled as a localized 
body force in the vicinity of the leading edge. We consider the output to be the entire velocity 
field: this is important for capturing the flow physics, and is often needed to represent cost functions 
used in optimal control design. Since the output is very high-dimensional, in ERA and balanced 
POD reduction procedures we use output projection described at the end of Section [2| projecting 
the velocity field onto the leading POD modes of the velocity snapshots obtained from the impulse 
response simulation. 



4.3 Reduced-order models 



ERA is applied to the full-dimensional linearized system to construct a reduced-order model. With a 
sampling period of 50 time steps, 400 adjacent pairs of Markov parameters, as in ( fT0| ), are collected 
from an impulse response simulation. Since these parameters are a projection of the velocity fields 
onto the leading POD modes, for an output projection of order m, the number of inner products 



required is 4m x 10 for construction of each H and H' (see Section 2.4 1 



For comparison, balanced POD is also used to compute the same reduced-order models. Adjoint 
simulations are performed with the POD modes as initial conditions to compute the matrix Y of (|5]>. 



13 



Steps in computing Approximate time (CPU hours) 



reduced-order models 


balanced POD 


ERA 


1 . Linearized impulse response 


2 


4 


2. Computation of POD modes 


2 


2 


3. Adjoint impulse responses 


30 




(10 in number) 






4. Computation of the Hankel matrix 


7 


0.2 


5. Singular value decomposition 


0.05 


0.05 


6. Computation of modes 


1 




7. Computation of models 


0.02 


0.02 



Table 1 : Comparison of the computational times required for various steps of the algorithms using 
balanced POD and ERA. The times are given for a 10-mode output projected system. The Hankel 
matrix is constructed using (a) 200 state-snapshots from each linearized and adjoint simulations for 
balanced POD, and (b) 400 Markov parameters (outputs) for ERA. 

The matrices X and Y are assembled by stacking 200 snapshots from the linearized and each of 
the adjoint simulations, and in turn, the generalized Hankel matrix H = Y*X is computed. For an 
output projection of order m, the number of inner products required to compute H is Am x 10 4 , 
which is 50 times more than that to compute H and H' in total for ERA. 

We also compare reduced-order models using standard POD modes, and ERA with pseudo-adjoint 
modes, as described in Section[3] The first 100 primal modes are used to compute the pseudo-adjoint 
modes.. 

For the given case, a comparison between the computational cost using ERA and using balanced 
POD is shown in Table [T] Results verify that ERA substantially improves computational efficiency 
in forming reduced-order models. 

Next, we compare the reduced-order models. Figure [3] shows the leading two primal modes and 
true adjoint modes from balanced POD, compared with the leading two pseudo-adjoint modes. The 
pseudo-adjoint modes look quite different from the true adjoint modes, and the flow structures ac- 
tually more closely resemble the leading primal modes. This result is not surprising, since the 
pseudo-adjoint modes are always linear combinations of the snapshots from the primal simulations, 
while the true adjoint modes are linear combinations of snapshots from adjoint simulations. Follow- 
ing the discussion in the last section, the poor approximation of the adjoint modes suggests that the 
pseudo-adjoint modes may produce poor reduced order models for this example, as we will verify 
below. 

Figure [4] shows the diagonal values of the controllability and observability Gramians, as well as the 
empirical Hankel singular values, for reduced-order models obtained from three different methods: 
ERA, balanced POD, and ERA with pseudo-adjoint modes. The models obtained using ERA are 
more accurate in the sense that the three sets of curves are almost indistinguishable, for all orders of 



14 







1 



1 

0.5 
a 
-0.5 
-1 



Balanced POD, mode 1 




-1 



1 



1 

0.5 

-0.5 
-1 



Balanced POD, adjoint mode 1 



1 




0.5 




a. o 




-0.5 




-1 






Pseudo- adjoint mode 1 



1 



Balanced POD. mode 2 




-1 







1 



Balanced POD, adjoint mode 2 

m 




1 



Pseudo- adjoint mode 2 



Figure 3: The first two primal and adjoint modes computed by using balanced POD, and the first 
two pseudo-adjoint modes computed by using ( fl~4| ) and the first 100 primal modes. Modes are 
illustrated using contour plots of the vorticity field. 
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Figure 4: Comparison of Gramians computed using (a) ERA, (b) balanced POD, and (c) ERA with 

pseudo-adjoint modes: The empirical Hankel singular values ( ) and the diagonal elements of 

the controllability ( , o) and observability ( , x ) Gramians with different order of modes 

(e.g., 4, 10, 20) in output projection. 
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output-projection. However, for balanced POD, the diagonal values of the observability Gramians 
are accurate only for certain leading modes, the number of which depends on and increases with 
the order of output projection. This inaccuracy can be attributed to a slight inaccuracy in the ad- 
joint formulation, which in turn results from an approximation in the multi-domain approach used 
to treat f arfield boundary conditions in the immersed boundary method of ; see J2 for more de- 
tails. Thus, ERA is advantageous as it does not need any adjoint simulations and results in more 
balanced Gramians. On the other hand, ERA with pseudo-adjoint modes generates poorly balanced 
controllability and observability Gramians, as shown in Figure ^c). This is because the leading 
primal modes and adjoint modes are supported very differently in the spatial domain, and thus the 
pseudo-adjoint modes, based on linear combination of leading primal modes, poorly approximate 
the true adjoint modes. 

4.4 Model performance 

We can quantify the performance of the various reduced-order models by computing error norms. 
One such measure is the 2-norm of the error between the impulse response of the full linearized 
system, denoted G(t), and that of a reduced order model with r modes, denoted by G r (t). We first 
compute the 2-norm of the error between the full system (with the entire velocity field as output) 
and the output-projected system of order 20, shown as the horizontal dashed line in Figure [5] This 
is the lower error bound for any reduced order model of the given output-projected system. Results 
shown in Figure [5] indicate that the first several low-order models obtained by ERA and balanced 
POD generate slightly different 2-norms of error, presumably because of the slight inaccuracy in 
the adjoint, mentioned previously. For most orders, however, they agree, and both error norms 
converge to the lower bound as the order of the model increases. By running more simulation tests, 
we observe that with higher-order output projections, ERA and balanced POD error norms converge 
to each other faster when the order of model increases. 

Figure [5] also shows the 2-norm error plots for models by ERA with pseudo-adjoint modes, using 
20-mode output projection, and for models computed using standard POD. Errors of models by 
ERA with pseudo-adjoint modes converge to the lower bound much slower than ERA/balanced 
POD. Errors of models by POD do not start converging until more than nearly 20 modes are used, 
and they converge to a larger error bound than ERA/balanced POD, again because POD models do 
not capture the input-output dynamics as well as balanced truncation based models. 

In the time domain, a comparison of the transient response to an impulsive disturbance is shown 
in Figure |6| in which the first output of the reduced-order model is plotted, for a 16-mode model 
determined by ERA, and for 30-mode models by POD and ERA with pseudo-adjoint modes . The 
16-mode ERA model already accurately predicts the response for all times. The higher-dimensional, 
30-mode models using POD and pseudo-adjoint modes are both stable, and perform reasonably 
well; however, they over-predict the response, particularly after time t « 80. 
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Figure 5: H2— norm of the error with increasing order of the reduced-order models: exact out- 
put of the output-projected system ( ); models obtained using balanced POD ( , □), 

ERA( , x ), ERA with pseudo-adjoint modes( — , V) , and POD( , A). A 20-mode output 

projection is used in ERA, balanced POD, and ERA with pseudo-adjoint modes. 




Figure 6: The first output, output a\, from the impulse-response simulation: results of full- 
simulation ( , o), compared with those of 16-mode reduced order model by ERA( , x), 30- 

mode model by ERA with pseudo-adjoint modes( — , V) , and 30-mode model by POD( , A) 

. A 20-mode output projection is used in ERA and ERA with pseudo-adjoint modes. 



18 



-50 



10 



10 u 10' 

frequency (Hz) 



10 



Figure 7: Singular- value plots: The full system ( ) and 30-mode models obtained using bal- 
anced POD ( ), ERA ( ), ERA with pseudo-adjoint modes( ), and POD( ), all 

with a 20-mode output projection. ERA and balanced POD models generate almost identical plots. 



We also compare the frequency response of reduced-order models to that of the full system, or more 
precisely, the full output-projected system. One way to represent the response of a single-input 
multiple-output system is by a singular-value plot, a plot of the maximum singular value of the 
transfer function matrix as a function of frequency. To generate this plot, a very long simulation 
of 5 x 10 5 time steps for the full system is performed, with a random input sampled from a uni- 
form distribution in the range (—0.5, 0.5). The output snapshots are projected onto leading POD 
modes. The magnitude of the transfer function is then computed from the cross spectrum of the 
input and outputs (using the Matlab command tf estimate). Finally, singular- value plots for the 
full output-projected systems are obtained, with a typical case shown as the red line in Figure [7] 

A typical set of singular- value plots of different reduced order models are presented in Figure [7] 
Results shown in the figure indicate that ERA and balanced POD 30-mode models, are almost 
identical, and are close to the corresponding full output-projected system. In comparison, Figure [7] 
also shows sigma plots for 30-mode models by ERA with pseudo-adjoint modes and by POD. Note 
that for computational feasibility, here the output of the POD model is the first twenty reduced 
states, i.e., the full-dimensional output of the POD model are projected onto the leading twenty 
POD modes. The frequency responses of the models by POD and ERA with pseudo-adjoint modes 
capture the resonant peak, but do not match well for frequencies far away from the resonant peak. 
These two models both generate spurious peaks in the frequency range of [0.1, 2]. 



5 Discussion 



We report that, theoretically, the eigensystem realization algorithm (ERA) and snapshot-based ap- 
proximate balanced truncation (balanced POD) produce exactly the same reduced order models. 
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This equivalence implies that ERA balances a pair of approximate Gramians and thus can be re- 
garded as an approximate balanced truncation method. Compared to balanced POD, the main 
features of ERA are that it does not require data from adjoint systems and therefore can be used 
with experimental data; furthermore, its construction of the generalized Hankel matrix is compu- 
tationally an order-of-magnitude cheaper than balanced POD. Numerical results indicate that ERA 
can be more accurate than balanced POD in practice, since there can be slight inaccuracies in the 
adjoint operator used with balanced POD. Balanced POD does have its own advantages, however: 
unlike ERA, it produces sets of bi-orthogonal modes that are useful for other purposes. Nonlinear 
models may be obtained by projection onto these modes; and parameters such as Reynolds number 
can be retained in the reduced-order models generated using these modes. Balanced POD has also 
been generalized for unstable systems. We also examine a modified ERA approach in which one 
constructs sets of bi-orthogonal modes without using adjoint information, using a matrix pseudo- 
inverse, as in [22]. Although this approach provides sets of bi-orthogonal modes (primal/pseudo- 
adjoint modes), in general it can not be regarded as an approximate balanced truncation method, 
since it does not balance the approximate Gramians. 

We have demonstrated the methods on an model problem consisting of a disturbance interacting with 
the flow past an inclined flat plate. As expected, balanced POD models perform nearly identically 
to ERA models. The small differences result because the adjoint simulation required for balanced 
POD is not a perfect adjoint at the discrete level. Both procedures work significantly better than 
standard POD models, or ERA models using pseudo-adjoint modes for projection. 

Finally, we emphasize that throughout, we have considered only stable, linear models. Possible fu- 
ture directions of this work include a generalization to unstable systems, and ultimately to nonlinear 
systems. 
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